fig3_Polarized_and_Naive

#polar 1

DE RNA

# month 3 comparisons 

Meth_vs_Nal_3<-FindMarkers(results_polar1, "Methadone_3","Naltrexone_3", min.pct = -Inf, min.diff.pct = -Inf, logfc.threshold = -Inf)
Meth_vs_Bup.Nalo_3<-FindMarkers(results_polar1, "Methadone_3","Bup.Nalo_3", min.pct = -Inf, min.diff.pct = -Inf, logfc.threshold = -Inf)
Bup.Nalo_vs_Nal_3<-FindMarkers(results_polar1, "Bup.Nalo_3","Naltrexone_3", min.pct = -Inf, min.diff.pct = -Inf, logfc.threshold = -Inf)

# month 0 comparisons 
Meth_vs_Nal_0<-FindMarkers(results_polar1, "Methadone_0","Naltrexone_0", min.pct = -Inf, min.diff.pct = -Inf, logfc.threshold = -Inf)
Meth_vs_Bup.Nalo_0<-FindMarkers(results_polar1, "Methadone_0","Bup.Nalo_0", min.pct = -Inf, min.diff.pct = -Inf, logfc.threshold = -Inf)
Bup.Nalo_vs_Nal_0<-FindMarkers(results_polar1, "Bup.Nalo_0","Naltrexone_0", min.pct = -Inf, min.diff.pct = -Inf, logfc.threshold = -Inf)

todo<-list(Meth_vs_Nal_3,Meth_vs_Bup.Nalo_3,Bup.Nalo_vs_Nal_3,Meth_vs_Nal_0,Meth_vs_Bup.Nalo_0,Bup.Nalo_vs_Nal_0 )
names(todo)<-c("Meth_vs_Nal_3","Meth_vs_Bup.Nalo_3","Bup.Nalo_vs_Nal_3","Meth_vs_Nal_0","Meth_vs_Bup.Nalo_0","Bup.Nalo_vs_Nal_0")

for(i in 1:length(todo)){
  todo[[i]]$p_val_adj<-p.adjust( todo[[i]]$p_val, "BH")
  print(VolPlot( todo[[i]], Title = names(todo)[[i]]))
}

## Warning: ggrepel: 16 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

## Warning: ggrepel: 69 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

## Warning: ggrepel: 8 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

## Warning: ggrepel: 63 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

## Warning: ggrepel: 8 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

pdf("~/gibbs/DOGMAMORPH/Ranalysis/Scripts/Figure Notebooks/rawFigs/fig3/A_F.pdf", width = 16, height = 9)
for(i in 1:length(todo)){
  print(VolPlot(todo[[i]], Title = names(todo)[[i]]))
}

## Warning: ggrepel: 16 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

## Warning: ggrepel: 69 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

## Warning: ggrepel: 9 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

## Warning: ggrepel: 63 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

## Warning: ggrepel: 8 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

dev.off()->.

#subset to just DE genes for these tables to avoid them being too unwieldy
DT::datatable(rownames=TRUE, filter="top", class='cell-border stripe', extensions = 'Buttons', options = list(dom = 'Bfrtip', buttons = c('copy', 'csv', 'excel', 'pdf', 'print')), data=subset(Meth_vs_Nal_3, p_val_adj<0.01))

DT::datatable(rownames=TRUE, filter="top", class='cell-border stripe', extensions = 'Buttons', options = list(dom = 'Bfrtip', buttons = c('copy', 'csv', 'excel', 'pdf', 'print')), data=subset(Meth_vs_Bup.Nalo_3, p_val_adj<0.01))

DT::datatable(rownames=TRUE, filter="top", class='cell-border stripe', extensions = 'Buttons', options = list(dom = 'Bfrtip', buttons = c('copy', 'csv', 'excel', 'pdf', 'print')), data=subset(Bup.Nalo_vs_Nal_3, p_val_adj<0.01))

DT::datatable(rownames=TRUE, filter="top", class='cell-border stripe', extensions = 'Buttons', options = list(dom = 'Bfrtip', buttons = c('copy', 'csv', 'excel', 'pdf', 'print')), data=subset(Meth_vs_Nal_0, p_val_adj<0.01))

DT::datatable(rownames=TRUE, filter="top", class='cell-border stripe', extensions = 'Buttons', options = list(dom = 'Bfrtip', buttons = c('copy', 'csv', 'excel', 'pdf', 'print')), data=subset(Meth_vs_Bup.Nalo_0, p_val_adj<0.01))

DT::datatable(rownames=TRUE, filter="top", class='cell-border stripe', extensions = 'Buttons', options = list(dom = 'Bfrtip', buttons = c('copy', 'csv', 'excel', 'pdf', 'print')), data=subset(Bup.Nalo_vs_Nal_0, p_val_adj<0.01))

GSEA

#grabbing hallmark as well as the curated, immune 
m_df_H<- msigdbr(species = "Homo sapiens", category = "H")
m_df_H<- rbind(msigdbr(species = "Homo sapiens", category = "C2"), m_df_H)
m_df_H<- rbind(msigdbr(species = "Homo sapiens", category = "C7"), m_df_H)
fgsea_sets<- m_df_H %>% split(x = .$gene_symbol, f = .$gs_name)

# month 3 comparisons 
Meth_vs_Nal_3<-FindMarkers(results_polar1, "Methadone_3","Naltrexone_3", min.pct = -Inf, min.diff.pct = -Inf, logfc.threshold = -Inf)
Bup.Nalo_vs_Meth_3<-FindMarkers(results_polar1, "Bup.Nalo_3","Methadone_3", min.pct = -Inf, min.diff.pct = -Inf, logfc.threshold = -Inf)
Bup.Nalo_vs_Nal_3<-FindMarkers(results_polar1, "Bup.Nalo_3","Naltrexone_3", min.pct = -Inf, min.diff.pct = -Inf, logfc.threshold = -Inf)

# month 0 comparisons 
Meth_vs_Nal_0<-FindMarkers(results_polar1, "Methadone_0","Naltrexone_0", min.pct = -Inf, min.diff.pct = -Inf, logfc.threshold = -Inf)
Bup.Nalo_vs_Meth_0<-FindMarkers(results_polar1, "Bup.Nalo_0","Methadone_0", min.pct = -Inf, min.diff.pct = -Inf, logfc.threshold = -Inf)
Bup.Nalo_vs_Nal_0<-FindMarkers(results_polar1, "Bup.Nalo_0","Naltrexone_0", min.pct = -Inf, min.diff.pct = -Inf, logfc.threshold = -Inf)

todo<-list(Meth_vs_Nal_3,Bup.Nalo_vs_Meth_3,Bup.Nalo_vs_Nal_3,Meth_vs_Nal_0,Bup.Nalo_vs_Meth_0,Bup.Nalo_vs_Nal_0 )
names(todo)<-c("Meth_vs_Nal_3","Bup.Nalo_vs_Meth_3","Bup.Nalo_vs_Nal_3","Meth_vs_Nal_0","Bup.Nalo_vs_Meth_0","Bup.Nalo_vs_Nal_0")

GSEAres<-list()
for (i in 1:length(todo)){
GSEAres[[i]]<-GSEA(todo[[i]], genesets = fgsea_sets)
GSEAres[[i]]<-GSEATable(GSEAwrap_out =GSEAres[[i]], gmt = fgsea_sets, name = names(todo)[[i]] )
}

## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, : There are ties in the preranked stats (7.31% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.

## [1] "start ranking"
## [1] "done ranking"

## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, : There are ties in the preranked stats (6.48% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.

## Warning in fgseaMultilevel(pathways = pathways, stats = stats, minSize =
## minSize, : There were 8 pathways for which P-values were not calculated
## properly due to unbalanced (positive and negative) gene-level statistic values.
## For such pathways pval, padj, NES, log2err are set to NA. You can try to
## increase the value of the argument nPermSimple (for example set it nPermSimple
## = 10000)

## [1] "start ranking"
## [1] "done ranking"

## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, : There are ties in the preranked stats (6.51% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.

## Warning in fgseaMultilevel(pathways = pathways, stats = stats, minSize =
## minSize, : There were 9 pathways for which P-values were not calculated
## properly due to unbalanced (positive and negative) gene-level statistic values.
## For such pathways pval, padj, NES, log2err are set to NA. You can try to
## increase the value of the argument nPermSimple (for example set it nPermSimple
## = 10000)

## [1] "start ranking"
## [1] "done ranking"

## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, : There are ties in the preranked stats (6.11% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.

## Warning in fgseaMultilevel(pathways = pathways, stats = stats, minSize =
## minSize, : There were 7 pathways for which P-values were not calculated
## properly due to unbalanced (positive and negative) gene-level statistic values.
## For such pathways pval, padj, NES, log2err are set to NA. You can try to
## increase the value of the argument nPermSimple (for example set it nPermSimple
## = 10000)

## Warning in fgseaMultilevel(pathways = pathways, stats = stats, minSize =
## minSize, : For some of the pathways the P-values were likely overestimated. For
## such pathways log2err is set to NA.

## [1] "start ranking"
## [1] "done ranking"

## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, : There are ties in the preranked stats (5.73% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.

## [1] "start ranking"
## [1] "done ranking"

## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, : There are ties in the preranked stats (4.5% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.

## Warning in fgseaMultilevel(pathways = pathways, stats = stats, minSize =
## minSize, : There were 16 pathways for which P-values were not calculated
## properly due to unbalanced (positive and negative) gene-level statistic values.
## For such pathways pval, padj, NES, log2err are set to NA. You can try to
## increase the value of the argument nPermSimple (for example set it nPermSimple
## = 10000)

## [1] "start ranking"
## [1] "done ranking"

GSEAres<-GSEAbig(listofGSEAtables = GSEAres)


to_plot<-c("GSE11057_CD4_EFF_MEM_VS_PBMC_UP","GSE8685_IL2_STARVED_VS_IL2_ACT_IL2_STARVED_CD4_TCELL_DN","BROWNE_INTERFERON_RESPONSIVE_GENES")
DT::datatable(rownames=TRUE, filter="top", class='cell-border stripe', extensions = 'Buttons', options = list(dom = 'Bfrtip', buttons = c('copy', 'csv', 'excel', 'pdf', 'print')), data=subset(GSEAres, padj<0.001))

## Warning in instance$preRenderHook(instance): It seems your data is too big for
## client-side DataTables. You may consider server-side processing:
## https://rstudio.github.io/DT/server.html

GSEAEnrichmentPlotComparison(to_plot[1], GSEAres, "BOTH", cols = IsleofDogs1)

GSEAEnrichmentPlotComparison(to_plot[2], GSEAres, "BOTH", cols = IsleofDogs1)

GSEAEnrichmentPlotComparison(to_plot[3], GSEAres, "BOTH", cols = IsleofDogs1)

pdf("~/gibbs/DOGMAMORPH/Ranalysis/Scripts/Figure Notebooks/rawFigs/fig3/G_I.pdf", width = 16, height = 9)
GSEAEnrichmentPlotComparison(to_plot[1], GSEAres, "BOTH", cols = IsleofDogs1)
GSEAEnrichmentPlotComparison(to_plot[2], GSEAres, "BOTH", cols = IsleofDogs1)
GSEAEnrichmentPlotComparison(to_plot[3], GSEAres, "BOTH", cols = IsleofDogs1)
dev.off()

## png 
##   2

DE CHROMVar

DefaultAssay(results_polar1)<-"chromvar"
# month 3 comparisons 

Meth_vs_Nal_3<-FindMarkers(results_polar1, "Methadone_3","Naltrexone_3", mean.fxn = rowMeans)
Meth_vs_Bup.Nalo_3<-FindMarkers(results_polar1, "Methadone_3","Bup.Nalo_3", mean.fxn = rowMeans)
Bup.Nalo_vs_Nal_3<-FindMarkers(results_polar1, "Bup.Nalo_3","Naltrexone_3", mean.fxn = rowMeans)

# month 0 comparisons 
Meth_vs_Nal_0<-FindMarkers(results_polar1, "Methadone_0","Naltrexone_0", mean.fxn = rowMeans)
Meth_vs_Bup.Nalo_0<-FindMarkers(results_polar1, "Methadone_0","Bup.Nalo_0", mean.fxn = rowMeans)
Bup.Nalo_vs_Nal_0<-FindMarkers(results_polar1, "Bup.Nalo_0","Naltrexone_0", mean.fxn = rowMeans)

Meth_vs_Nal_3<-FixMotifID(Meth_vs_Nal_3, results_polar1)
Meth_vs_Bup.Nalo_3<-FixMotifID(Meth_vs_Bup.Nalo_3, results_polar1)
Bup.Nalo_vs_Nal_3<-FixMotifID(Bup.Nalo_vs_Nal_3, results_polar1)
Meth_vs_Nal_0<-FixMotifID(Meth_vs_Nal_0, results_polar1)
Meth_vs_Bup.Nalo_0<-FixMotifID(Meth_vs_Bup.Nalo_0, results_polar1)
Bup.Nalo_vs_Nal_0<-FixMotifID(Bup.Nalo_vs_Nal_0, results_polar1)

DT::datatable(rownames=TRUE, filter="top", class='cell-border stripe', extensions = 'Buttons', options = list(dom = 'Bfrtip', buttons = c('copy', 'csv', 'excel', 'pdf', 'print')), data=Meth_vs_Nal_3)

DT::datatable(rownames=TRUE, filter="top", class='cell-border stripe', extensions = 'Buttons', options = list(dom = 'Bfrtip', buttons = c('copy', 'csv', 'excel', 'pdf', 'print')), data=Meth_vs_Bup.Nalo_3)

DT::datatable(rownames=TRUE, filter="top", class='cell-border stripe', extensions = 'Buttons', options = list(dom = 'Bfrtip', buttons = c('copy', 'csv', 'excel', 'pdf', 'print')), data=Bup.Nalo_vs_Nal_3)

DT::datatable(rownames=TRUE, filter="top", class='cell-border stripe', extensions = 'Buttons', options = list(dom = 'Bfrtip', buttons = c('copy', 'csv', 'excel', 'pdf', 'print')), data=Meth_vs_Nal_0)

DT::datatable(rownames=TRUE, filter="top", class='cell-border stripe', extensions = 'Buttons', options = list(dom = 'Bfrtip', buttons = c('copy', 'csv', 'excel', 'pdf', 'print')), data=Meth_vs_Bup.Nalo_0)

DT::datatable(rownames=TRUE, filter="top", class='cell-border stripe', extensions = 'Buttons', options = list(dom = 'Bfrtip', buttons = c('copy', 'csv', 'excel', 'pdf', 'print')), data=Bup.Nalo_vs_Nal_0)

Idents(results_polar1)<-factor(Idents(results_polar1),levels = c("Methadone_0","Bup.Nalo_0","Naltrexone_0","Methadone_3", "Bup.Nalo_3", "Naltrexone_3"))

BottleRocket3<-c("Methadone_0" = "#3b4357", "Bup.Nalo_0" = "#cb2314","Naltrexone_0"= "#fad510",
                 "Methadone_3" = "#273046", "Bup.Nalo_3" = "#792a2d","Naltrexone_3"= "#e37c12")

pdf("~/gibbs/DOGMAMORPH/Ranalysis/Scripts/Figure Notebooks/rawFigs/fig3/J_K.pdf", width = 16, height = 9)
VlnPlot(results_polar1, "MA1143.1", pt.size = 0.1 )+ggtitle("FOSL1::JUND")+scale_fill_manual(values=BottleRocket3)
VlnPlot(results_polar1, "MA0605.2", pt.size = 0.1 )+ggtitle("ATF3")+scale_fill_manual(values=BottleRocket3)
dev.off()

## png 
##   2

polar 2

results_polar2<-subset(results, merged_clusters=="MemoryT_Polarized_2")
Idents(results_polar2)<-results_polar2$T_Tp
DefaultAssay(results_polar2)<-"RNA"

DE RNA

# month 3 comparisons 

Meth_vs_Nal_3<-FindMarkers(results_polar2, "Methadone_3","Naltrexone_3", min.pct = -Inf, min.diff.pct = -Inf, logfc.threshold = -Inf)
Meth_vs_Bup.Nalo_3<-FindMarkers(results_polar2, "Methadone_3","Bup.Nalo_3", min.pct = -Inf, min.diff.pct = -Inf, logfc.threshold = -Inf)
Bup.Nalo_vs_Nal_3<-FindMarkers(results_polar2, "Bup.Nalo_3","Naltrexone_3", min.pct = -Inf, min.diff.pct = -Inf, logfc.threshold = -Inf)

# month 0 comparisons 
Meth_vs_Nal_0<-FindMarkers(results_polar2, "Methadone_0","Naltrexone_0", min.pct = -Inf, min.diff.pct = -Inf, logfc.threshold = -Inf)
Meth_vs_Bup.Nalo_0<-FindMarkers(results_polar2, "Methadone_0","Bup.Nalo_0", min.pct = -Inf, min.diff.pct = -Inf, logfc.threshold = -Inf)
Bup.Nalo_vs_Nal_0<-FindMarkers(results_polar2, "Bup.Nalo_0","Naltrexone_0", min.pct = -Inf, min.diff.pct = -Inf, logfc.threshold = -Inf)

todo<-list(Meth_vs_Nal_3,Meth_vs_Bup.Nalo_3,Bup.Nalo_vs_Nal_3,Meth_vs_Nal_0,Meth_vs_Bup.Nalo_0,Bup.Nalo_vs_Nal_0 )
names(todo)<-c("Meth_vs_Nal_3","Meth_vs_Bup.Nalo_3","Bup.Nalo_vs_Nal_3","Meth_vs_Nal_0","Meth_vs_Bup.Nalo_0","Bup.Nalo_vs_Nal_0")

for(i in 1:length(todo)){
  todo[[i]]$p_val_adj<-p.adjust( todo[[i]]$p_val, "BH")
  print(VolPlot( todo[[i]], Title = names(todo)[[i]]))
}

## Warning: ggrepel: 42 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

## Warning: ggrepel: 26 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

## Warning: ggrepel: 1 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

pdf("~/gibbs/DOGMAMORPH/Ranalysis/Scripts/Figure Notebooks/rawFigs/fig3/L_Q.pdf", width = 16, height = 9)
for(i in 1:length(todo)){
  print(VolPlot(todo[[i]], Title = names(todo)[[i]]))
}

## Warning: ggrepel: 41 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

## Warning: ggrepel: 25 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

## Warning: ggrepel: 1 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

dev.off()->.

#subset to just DE genes for these tables to avoid them being too unwieldy
DT::datatable(rownames=TRUE, filter="top", class='cell-border stripe', extensions = 'Buttons', options = list(dom = 'Bfrtip', buttons = c('copy', 'csv', 'excel', 'pdf', 'print')), data=subset(Meth_vs_Nal_3, p_val_adj<0.01))

DT::datatable(rownames=TRUE, filter="top", class='cell-border stripe', extensions = 'Buttons', options = list(dom = 'Bfrtip', buttons = c('copy', 'csv', 'excel', 'pdf', 'print')), data=subset(Meth_vs_Bup.Nalo_3, p_val_adj<0.01))

DT::datatable(rownames=TRUE, filter="top", class='cell-border stripe', extensions = 'Buttons', options = list(dom = 'Bfrtip', buttons = c('copy', 'csv', 'excel', 'pdf', 'print')), data=subset(Bup.Nalo_vs_Nal_3, p_val_adj<0.01))

DT::datatable(rownames=TRUE, filter="top", class='cell-border stripe', extensions = 'Buttons', options = list(dom = 'Bfrtip', buttons = c('copy', 'csv', 'excel', 'pdf', 'print')), data=subset(Meth_vs_Nal_0, p_val_adj<0.01))

DT::datatable(rownames=TRUE, filter="top", class='cell-border stripe', extensions = 'Buttons', options = list(dom = 'Bfrtip', buttons = c('copy', 'csv', 'excel', 'pdf', 'print')), data=subset(Meth_vs_Bup.Nalo_0, p_val_adj<0.01))

DT::datatable(rownames=TRUE, filter="top", class='cell-border stripe', extensions = 'Buttons', options = list(dom = 'Bfrtip', buttons = c('copy', 'csv', 'excel', 'pdf', 'print')), data=subset(Bup.Nalo_vs_Nal_0, p_val_adj<0.01))

GSEA

#grabbing hallmark as well as the curated, immune 
m_df_H<- msigdbr(species = "Homo sapiens", category = "H")
m_df_H<- rbind(msigdbr(species = "Homo sapiens", category = "C2"), m_df_H)
m_df_H<- rbind(msigdbr(species = "Homo sapiens", category = "C7"), m_df_H)
fgsea_sets<- m_df_H %>% split(x = .$gene_symbol, f = .$gs_name)

# month 3 comparisons 
Meth_vs_Nal_3<-FindMarkers(results_polar2, "Methadone_3","Naltrexone_3", min.pct = -Inf, min.diff.pct = -Inf, logfc.threshold = -Inf)
Bup.Nalo_vs_Meth_3<-FindMarkers(results_polar2, "Bup.Nalo_3","Methadone_3", min.pct = -Inf, min.diff.pct = -Inf, logfc.threshold = -Inf)
Bup.Nalo_vs_Nal_3<-FindMarkers(results_polar2, "Bup.Nalo_3","Naltrexone_3", min.pct = -Inf, min.diff.pct = -Inf, logfc.threshold = -Inf)

# month 0 comparisons 
Meth_vs_Nal_0<-FindMarkers(results_polar2, "Methadone_0","Naltrexone_0", min.pct = -Inf, min.diff.pct = -Inf, logfc.threshold = -Inf)
Bup.Nalo_vs_Meth_0<-FindMarkers(results_polar2, "Bup.Nalo_0","Methadone_0", min.pct = -Inf, min.diff.pct = -Inf, logfc.threshold = -Inf)
Bup.Nalo_vs_Nal_0<-FindMarkers(results_polar2, "Bup.Nalo_0","Naltrexone_0", min.pct = -Inf, min.diff.pct = -Inf, logfc.threshold = -Inf)

todo<-list(Meth_vs_Nal_3,Bup.Nalo_vs_Meth_3,Bup.Nalo_vs_Nal_3,Meth_vs_Nal_0,Bup.Nalo_vs_Meth_0,Bup.Nalo_vs_Nal_0 )
names(todo)<-c("Meth_vs_Nal_3","Bup.Nalo_vs_Meth_3","Bup.Nalo_vs_Nal_3","Meth_vs_Nal_0","Bup.Nalo_vs_Meth_0","Bup.Nalo_vs_Nal_0")

GSEAres<-list()
for (i in 1:length(todo)){
GSEAres[[i]]<-GSEA(todo[[i]], genesets = fgsea_sets)
GSEAres[[i]]<-GSEATable(GSEAwrap_out =GSEAres[[i]], gmt = fgsea_sets, name = names(todo)[[i]] )
}

## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, : There are ties in the preranked stats (6.94% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.

## [1] "start ranking"
## [1] "done ranking"

## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, : There are ties in the preranked stats (6.26% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.

## [1] "start ranking"
## [1] "done ranking"

## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, : There are ties in the preranked stats (6.14% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.

## [1] "start ranking"
## [1] "done ranking"

## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, : There are ties in the preranked stats (6.59% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.

## Warning in fgseaMultilevel(pathways = pathways, stats = stats, minSize =
## minSize, : There were 7 pathways for which P-values were not calculated
## properly due to unbalanced (positive and negative) gene-level statistic values.
## For such pathways pval, padj, NES, log2err are set to NA. You can try to
## increase the value of the argument nPermSimple (for example set it nPermSimple
## = 10000)

## Warning in fgseaMultilevel(pathways = pathways, stats = stats, minSize =
## minSize, : For some of the pathways the P-values were likely overestimated. For
## such pathways log2err is set to NA.

## [1] "start ranking"
## [1] "done ranking"

## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, : There are ties in the preranked stats (5.98% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.

## Warning in fgseaMultilevel(pathways = pathways, stats = stats, minSize =
## minSize, : There were 4 pathways for which P-values were not calculated
## properly due to unbalanced (positive and negative) gene-level statistic values.
## For such pathways pval, padj, NES, log2err are set to NA. You can try to
## increase the value of the argument nPermSimple (for example set it nPermSimple
## = 10000)

## [1] "start ranking"
## [1] "done ranking"

## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, : There are ties in the preranked stats (4.84% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.

## [1] "start ranking"
## [1] "done ranking"

GSEAres<-GSEAbig(listofGSEAtables = GSEAres)


to_plot<-c("GSE21670_TGFB_VS_IL6_TREATED_STAT3_KO_CD4_TCELL_DN","BROCKE_APOPTOSIS_REVERSED_BY_IL6","PHONG_TNF_TARGETS_UP")
DT::datatable(rownames=TRUE, filter="top", class='cell-border stripe', extensions = 'Buttons', options = list(dom = 'Bfrtip', buttons = c('copy', 'csv', 'excel', 'pdf', 'print')), data=subset(GSEAres, padj<0.001))

## Warning in instance$preRenderHook(instance): It seems your data is too big for
## client-side DataTables. You may consider server-side processing:
## https://rstudio.github.io/DT/server.html

GSEAEnrichmentPlotComparison(to_plot[1], GSEAres, "BOTH", cols = IsleofDogs1)

GSEAEnrichmentPlotComparison(to_plot[2], GSEAres, "BOTH", cols = IsleofDogs1)

GSEAEnrichmentPlotComparison(to_plot[3], GSEAres, "BOTH", cols = IsleofDogs1)

pdf("~/gibbs/DOGMAMORPH/Ranalysis/Scripts/Figure Notebooks/rawFigs/fig3/R_T.pdf", width = 16, height = 9)
GSEAEnrichmentPlotComparison(to_plot[1], GSEAres, "BOTH", cols = IsleofDogs1)
GSEAEnrichmentPlotComparison(to_plot[2], GSEAres, "BOTH", cols = IsleofDogs1)
GSEAEnrichmentPlotComparison(to_plot[3], GSEAres, "BOTH", cols = IsleofDogs1)
dev.off()

## png 
##   2

DE CHROMVar

DefaultAssay(results_polar2)<-"chromvar"
# month 3 comparisons 

Meth_vs_Nal_3<-FindMarkers(results_polar2, "Methadone_3","Naltrexone_3", mean.fxn = rowMeans)
Meth_vs_Bup.Nalo_3<-FindMarkers(results_polar2, "Methadone_3","Bup.Nalo_3", mean.fxn = rowMeans)
Bup.Nalo_vs_Nal_3<-FindMarkers(results_polar2, "Bup.Nalo_3","Naltrexone_3", mean.fxn = rowMeans)

# month 0 comparisons 
Meth_vs_Nal_0<-FindMarkers(results_polar2, "Methadone_0","Naltrexone_0", mean.fxn = rowMeans)
Meth_vs_Bup.Nalo_0<-FindMarkers(results_polar2, "Methadone_0","Bup.Nalo_0", mean.fxn = rowMeans)
Bup.Nalo_vs_Nal_0<-FindMarkers(results_polar2, "Bup.Nalo_0","Naltrexone_0", mean.fxn = rowMeans)

Meth_vs_Nal_3<-FixMotifID(Meth_vs_Nal_3, results_polar2)
Meth_vs_Bup.Nalo_3<-FixMotifID(Meth_vs_Bup.Nalo_3, results_polar2)
Bup.Nalo_vs_Nal_3<-FixMotifID(Bup.Nalo_vs_Nal_3, results_polar2)
Meth_vs_Nal_0<-FixMotifID(Meth_vs_Nal_0, results_polar2)
Meth_vs_Bup.Nalo_0<-FixMotifID(Meth_vs_Bup.Nalo_0, results_polar2)
Bup.Nalo_vs_Nal_0<-FixMotifID(Bup.Nalo_vs_Nal_0, results_polar2)

DT::datatable(rownames=TRUE, filter="top", class='cell-border stripe', extensions = 'Buttons', options = list(dom = 'Bfrtip', buttons = c('copy', 'csv', 'excel', 'pdf', 'print')), data=Meth_vs_Nal_3)

DT::datatable(rownames=TRUE, filter="top", class='cell-border stripe', extensions = 'Buttons', options = list(dom = 'Bfrtip', buttons = c('copy', 'csv', 'excel', 'pdf', 'print')), data=Meth_vs_Bup.Nalo_3)

DT::datatable(rownames=TRUE, filter="top", class='cell-border stripe', extensions = 'Buttons', options = list(dom = 'Bfrtip', buttons = c('copy', 'csv', 'excel', 'pdf', 'print')), data=Bup.Nalo_vs_Nal_3)

DT::datatable(rownames=TRUE, filter="top", class='cell-border stripe', extensions = 'Buttons', options = list(dom = 'Bfrtip', buttons = c('copy', 'csv', 'excel', 'pdf', 'print')), data=Meth_vs_Nal_0)

DT::datatable(rownames=TRUE, filter="top", class='cell-border stripe', extensions = 'Buttons', options = list(dom = 'Bfrtip', buttons = c('copy', 'csv', 'excel', 'pdf', 'print')), data=Meth_vs_Bup.Nalo_0)

DT::datatable(rownames=TRUE, filter="top", class='cell-border stripe', extensions = 'Buttons', options = list(dom = 'Bfrtip', buttons = c('copy', 'csv', 'excel', 'pdf', 'print')), data=Bup.Nalo_vs_Nal_0)

Idents(results_polar2)<-factor(Idents(results_polar2),levels = c("Methadone_0","Bup.Nalo_0","Naltrexone_0","Methadone_3", "Bup.Nalo_3", "Naltrexone_3"))

BottleRocket3<-c("Methadone_0" = "#3b4357", "Bup.Nalo_0" = "#cb2314","Naltrexone_0"= "#fad510",
                 "Methadone_3" = "#273046", "Bup.Nalo_3" = "#792a2d","Naltrexone_3"= "#e37c12")

VlnPlot(results_polar2, "MA0688.1", pt.size = 0.1 )+ggtitle("TBX2")+scale_fill_manual(values=BottleRocket3)

VlnPlot(results_polar2, "MA0144.2", pt.size = 0.1 )+ggtitle("STAT3")+scale_fill_manual(values=BottleRocket3)

pdf("~/gibbs/DOGMAMORPH/Ranalysis/Scripts/Figure Notebooks/rawFigs/fig3/Q_R.pdf", width = 16, height = 9)
VlnPlot(results_polar2, "MA0688.1", pt.size = 0.1 )+ggtitle("TBX2")+scale_fill_manual(values=BottleRocket3)
VlnPlot(results_polar2, "MA0144.2", pt.size = 0.1 )+ggtitle("STAT3")+scale_fill_manual(values=BottleRocket3)
dev.off()

## png 
##   2

Naive cells

results_naive<-subset(results, merged_clusters=="NaiveT_1")
Idents(results_naive)<-results_naive$T_Tp
DefaultAssay(results_naive)<-"RNA"

DE RNA

# month 3 comparisons 

Meth_vs_Nal_3<-FindMarkers(results_naive, "Methadone_3","Naltrexone_3", min.pct = -Inf, min.diff.pct = -Inf, logfc.threshold = -Inf)
Meth_vs_Bup.Nalo_3<-FindMarkers(results_naive, "Methadone_3","Bup.Nalo_3", min.pct = -Inf, min.diff.pct = -Inf, logfc.threshold = -Inf)
Bup.Nalo_vs_Nal_3<-FindMarkers(results_naive, "Bup.Nalo_3","Naltrexone_3", min.pct = -Inf, min.diff.pct = -Inf, logfc.threshold = -Inf)

# month 0 comparisons 
Meth_vs_Nal_0<-FindMarkers(results_naive, "Methadone_0","Naltrexone_0", min.pct = -Inf, min.diff.pct = -Inf, logfc.threshold = -Inf)
Meth_vs_Bup.Nalo_0<-FindMarkers(results_naive, "Methadone_0","Bup.Nalo_0", min.pct = -Inf, min.diff.pct = -Inf, logfc.threshold = -Inf)
Bup.Nalo_vs_Nal_0<-FindMarkers(results_naive, "Bup.Nalo_0","Naltrexone_0", min.pct = -Inf, min.diff.pct = -Inf, logfc.threshold = -Inf)

todo<-list(Meth_vs_Nal_3,Meth_vs_Bup.Nalo_3,Bup.Nalo_vs_Nal_3,Meth_vs_Nal_0,Meth_vs_Bup.Nalo_0,Bup.Nalo_vs_Nal_0 )
names(todo)<-c("Meth_vs_Nal_3","Meth_vs_Bup.Nalo_3","Bup.Nalo_vs_Nal_3","Meth_vs_Nal_0","Meth_vs_Bup.Nalo_0","Bup.Nalo_vs_Nal_0")

for(i in 1:length(todo)){
  todo[[i]]$p_val_adj<-p.adjust( todo[[i]]$p_val, "BH")
  print(VolPlot( todo[[i]], Title = names(todo)[[i]]))
}

## Warning: ggrepel: 45 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

## Warning: ggrepel: 5 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

## Warning: ggrepel: 26 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

## Warning: ggrepel: 6 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

## Warning: ggrepel: 99 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

## Warning: ggrepel: 9 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

pdf("~/gibbs/DOGMAMORPH/Ranalysis/Scripts/Figure Notebooks/rawFigs/fig3/R_W.pdf", width = 16, height = 9)
for(i in 1:length(todo)){
  print(VolPlot(todo[[i]], Title = names(todo)[[i]]))
}

## Warning: ggrepel: 44 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

## Warning: ggrepel: 5 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

## Warning: ggrepel: 26 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

## Warning: ggrepel: 7 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

## Warning: ggrepel: 98 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

## Warning: ggrepel: 10 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

dev.off()->.

#subset to just DE genes for these tables to avoid them being too unwieldy
DT::datatable(rownames=TRUE, filter="top", class='cell-border stripe', extensions = 'Buttons', options = list(dom = 'Bfrtip', buttons = c('copy', 'csv', 'excel', 'pdf', 'print')), data=subset(Meth_vs_Nal_3, p_val_adj<0.01))

DT::datatable(rownames=TRUE, filter="top", class='cell-border stripe', extensions = 'Buttons', options = list(dom = 'Bfrtip', buttons = c('copy', 'csv', 'excel', 'pdf', 'print')), data=subset(Meth_vs_Bup.Nalo_3, p_val_adj<0.01))

DT::datatable(rownames=TRUE, filter="top", class='cell-border stripe', extensions = 'Buttons', options = list(dom = 'Bfrtip', buttons = c('copy', 'csv', 'excel', 'pdf', 'print')), data=subset(Bup.Nalo_vs_Nal_3, p_val_adj<0.01))

DT::datatable(rownames=TRUE, filter="top", class='cell-border stripe', extensions = 'Buttons', options = list(dom = 'Bfrtip', buttons = c('copy', 'csv', 'excel', 'pdf', 'print')), data=subset(Meth_vs_Nal_0, p_val_adj<0.01))

DT::datatable(rownames=TRUE, filter="top", class='cell-border stripe', extensions = 'Buttons', options = list(dom = 'Bfrtip', buttons = c('copy', 'csv', 'excel', 'pdf', 'print')), data=subset(Meth_vs_Bup.Nalo_0, p_val_adj<0.01))

DT::datatable(rownames=TRUE, filter="top", class='cell-border stripe', extensions = 'Buttons', options = list(dom = 'Bfrtip', buttons = c('copy', 'csv', 'excel', 'pdf', 'print')), data=subset(Bup.Nalo_vs_Nal_0, p_val_adj<0.01))

GSEA

#grabbing hallmark as well as the curated, immune 
m_df_H<- msigdbr(species = "Homo sapiens", category = "H")
m_df_H<- rbind(msigdbr(species = "Homo sapiens", category = "C2"), m_df_H)
m_df_H<- rbind(msigdbr(species = "Homo sapiens", category = "C7"), m_df_H)
fgsea_sets<- m_df_H %>% split(x = .$gene_symbol, f = .$gs_name)

# month 3 comparisons 
Meth_vs_Nal_3<-FindMarkers(results_naive, "Methadone_3","Naltrexone_3", min.pct = -Inf, min.diff.pct = -Inf, logfc.threshold = -Inf)
Bup.Nalo_vs_Meth_3<-FindMarkers(results_naive, "Bup.Nalo_3","Methadone_3", min.pct = -Inf, min.diff.pct = -Inf, logfc.threshold = -Inf)
Bup.Nalo_vs_Nal_3<-FindMarkers(results_naive, "Bup.Nalo_3","Naltrexone_3", min.pct = -Inf, min.diff.pct = -Inf, logfc.threshold = -Inf)

# month 0 comparisons 
Meth_vs_Nal_0<-FindMarkers(results_naive, "Methadone_0","Naltrexone_0", min.pct = -Inf, min.diff.pct = -Inf, logfc.threshold = -Inf)
Bup.Nalo_vs_Meth_0<-FindMarkers(results_naive, "Bup.Nalo_0","Methadone_0", min.pct = -Inf, min.diff.pct = -Inf, logfc.threshold = -Inf)
Bup.Nalo_vs_Nal_0<-FindMarkers(results_naive, "Bup.Nalo_0","Naltrexone_0", min.pct = -Inf, min.diff.pct = -Inf, logfc.threshold = -Inf)

todo<-list(Meth_vs_Nal_3,Bup.Nalo_vs_Meth_3,Bup.Nalo_vs_Nal_3,Meth_vs_Nal_0,Bup.Nalo_vs_Meth_0,Bup.Nalo_vs_Nal_0 )
names(todo)<-c("Meth_vs_Nal_3","Bup.Nalo_vs_Meth_3","Bup.Nalo_vs_Nal_3","Meth_vs_Nal_0","Bup.Nalo_vs_Meth_0","Bup.Nalo_vs_Nal_0")

GSEAres<-list()
for (i in 1:length(todo)){
GSEAres[[i]]<-GSEA(todo[[i]], genesets = fgsea_sets)
GSEAres[[i]]<-GSEATable(GSEAwrap_out =GSEAres[[i]], gmt = fgsea_sets, name = names(todo)[[i]] )
}

## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, : There are ties in the preranked stats (7.38% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.

## [1] "start ranking"
## [1] "done ranking"

## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, : There are ties in the preranked stats (6.38% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.

## Warning in fgseaMultilevel(pathways = pathways, stats = stats, minSize =
## minSize, : There were 2 pathways for which P-values were not calculated
## properly due to unbalanced (positive and negative) gene-level statistic values.
## For such pathways pval, padj, NES, log2err are set to NA. You can try to
## increase the value of the argument nPermSimple (for example set it nPermSimple
## = 10000)

## [1] "start ranking"
## [1] "done ranking"

## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, : There are ties in the preranked stats (6.75% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.

## [1] "start ranking"
## [1] "done ranking"

## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, : There are ties in the preranked stats (6.02% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.

## Warning in fgseaMultilevel(pathways = pathways, stats = stats, minSize =
## minSize, : There were 4 pathways for which P-values were not calculated
## properly due to unbalanced (positive and negative) gene-level statistic values.
## For such pathways pval, padj, NES, log2err are set to NA. You can try to
## increase the value of the argument nPermSimple (for example set it nPermSimple
## = 10000)

## [1] "start ranking"
## [1] "done ranking"

## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, : There are ties in the preranked stats (5.27% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.

## Warning in fgseaMultilevel(pathways = pathways, stats = stats, minSize =
## minSize, : There were 12 pathways for which P-values were not calculated
## properly due to unbalanced (positive and negative) gene-level statistic values.
## For such pathways pval, padj, NES, log2err are set to NA. You can try to
## increase the value of the argument nPermSimple (for example set it nPermSimple
## = 10000)

## [1] "start ranking"
## [1] "done ranking"

## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, : There are ties in the preranked stats (4.2% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.

## Warning in fgseaMultilevel(pathways = pathways, stats = stats, minSize =
## minSize, : There were 34 pathways for which P-values were not calculated
## properly due to unbalanced (positive and negative) gene-level statistic values.
## For such pathways pval, padj, NES, log2err are set to NA. You can try to
## increase the value of the argument nPermSimple (for example set it nPermSimple
## = 10000)

## [1] "start ranking"
## [1] "done ranking"

GSEAres<-GSEAbig(listofGSEAtables = GSEAres)


to_plot<-c("GSE2770_TGFB_AND_IL4_ACT_VS_ACT_CD4_TCELL_2H_UP","HALLMARK_INTERFERON_GAMMA_RESPONSE","PHONG_TNF_TARGETS_UP")
DT::datatable(rownames=TRUE, filter="top", class='cell-border stripe', extensions = 'Buttons', options = list(dom = 'Bfrtip', buttons = c('copy', 'csv', 'excel', 'pdf', 'print')), data=subset(GSEAres, padj<0.001))

## Warning in instance$preRenderHook(instance): It seems your data is too big for
## client-side DataTables. You may consider server-side processing:
## https://rstudio.github.io/DT/server.html

GSEAEnrichmentPlotComparison(to_plot[1], GSEAres, "BOTH", cols = IsleofDogs1)

GSEAEnrichmentPlotComparison(to_plot[2], GSEAres, "BOTH", cols = IsleofDogs1)

GSEAEnrichmentPlotComparison(to_plot[3], GSEAres, "BOTH", cols = IsleofDogs1)

pdf("~/gibbs/DOGMAMORPH/Ranalysis/Scripts/Figure Notebooks/rawFigs/fig3/X_Z.pdf", width = 16, height = 9)
GSEAEnrichmentPlotComparison(to_plot[1], GSEAres, "BOTH", cols = IsleofDogs1)
GSEAEnrichmentPlotComparison(to_plot[2], GSEAres, "BOTH", cols = IsleofDogs1)
GSEAEnrichmentPlotComparison(to_plot[3], GSEAres, "BOTH", cols = IsleofDogs1)
dev.off()

## png 
##   2

DE CHROMVar

DefaultAssay(results_naive)<-"chromvar"
# month 3 comparisons 

Meth_vs_Nal_3<-FindMarkers(results_naive, "Methadone_3","Naltrexone_3", mean.fxn = rowMeans)
Meth_vs_Bup.Nalo_3<-FindMarkers(results_naive, "Methadone_3","Bup.Nalo_3", mean.fxn = rowMeans)
Bup.Nalo_vs_Nal_3<-FindMarkers(results_naive, "Bup.Nalo_3","Naltrexone_3", mean.fxn = rowMeans)

# month 0 comparisons 
Meth_vs_Nal_0<-FindMarkers(results_naive, "Methadone_0","Naltrexone_0", mean.fxn = rowMeans)
Meth_vs_Bup.Nalo_0<-FindMarkers(results_naive, "Methadone_0","Bup.Nalo_0", mean.fxn = rowMeans)
Bup.Nalo_vs_Nal_0<-FindMarkers(results_naive, "Bup.Nalo_0","Naltrexone_0", mean.fxn = rowMeans)

Meth_vs_Nal_3<-FixMotifID(Meth_vs_Nal_3, results_naive)
Meth_vs_Bup.Nalo_3<-FixMotifID(Meth_vs_Bup.Nalo_3, results_naive)
Bup.Nalo_vs_Nal_3<-FixMotifID(Bup.Nalo_vs_Nal_3, results_naive)
Meth_vs_Nal_0<-FixMotifID(Meth_vs_Nal_0, results_naive)
Meth_vs_Bup.Nalo_0<-FixMotifID(Meth_vs_Bup.Nalo_0, results_naive)
Bup.Nalo_vs_Nal_0<-FixMotifID(Bup.Nalo_vs_Nal_0, results_naive)

DT::datatable(rownames=TRUE, filter="top", class='cell-border stripe', extensions = 'Buttons', options = list(dom = 'Bfrtip', buttons = c('copy', 'csv', 'excel', 'pdf', 'print')), data=Meth_vs_Nal_3)

DT::datatable(rownames=TRUE, filter="top", class='cell-border stripe', extensions = 'Buttons', options = list(dom = 'Bfrtip', buttons = c('copy', 'csv', 'excel', 'pdf', 'print')), data=Meth_vs_Bup.Nalo_3)

DT::datatable(rownames=TRUE, filter="top", class='cell-border stripe', extensions = 'Buttons', options = list(dom = 'Bfrtip', buttons = c('copy', 'csv', 'excel', 'pdf', 'print')), data=Bup.Nalo_vs_Nal_3)

DT::datatable(rownames=TRUE, filter="top", class='cell-border stripe', extensions = 'Buttons', options = list(dom = 'Bfrtip', buttons = c('copy', 'csv', 'excel', 'pdf', 'print')), data=Meth_vs_Nal_0)

DT::datatable(rownames=TRUE, filter="top", class='cell-border stripe', extensions = 'Buttons', options = list(dom = 'Bfrtip', buttons = c('copy', 'csv', 'excel', 'pdf', 'print')), data=Meth_vs_Bup.Nalo_0)

DT::datatable(rownames=TRUE, filter="top", class='cell-border stripe', extensions = 'Buttons', options = list(dom = 'Bfrtip', buttons = c('copy', 'csv', 'excel', 'pdf', 'print')), data=Bup.Nalo_vs_Nal_0)

Idents(results_naive)<-factor(Idents(results_naive),levels = c("Methadone_0","Bup.Nalo_0","Naltrexone_0","Methadone_3", "Bup.Nalo_3", "Naltrexone_3"))

BottleRocket3<-c("Methadone_0" = "#3b4357", "Bup.Nalo_0" = "#cb2314","Naltrexone_0"= "#fad510",
                 "Methadone_3" = "#273046", "Bup.Nalo_3" = "#792a2d","Naltrexone_3"= "#e37c12")

VlnPlot(results_naive, "MA0101.1", pt.size = 0.1 )+ggtitle("RELA")+scale_fill_manual(values=BottleRocket3)

VlnPlot(results_naive, "MA1143.1", pt.size = 0.1 )+ggtitle("FOSL1::JUND")+scale_fill_manual(values=BottleRocket3)

VlnPlot(results_naive, "MA1512.1", pt.size = 0.1 )+ggtitle("KLF11")+scale_fill_manual(values=BottleRocket3)

pdf("~/gibbs/DOGMAMORPH/Ranalysis/Scripts/Figure Notebooks/rawFigs/fig3/AA_AC.pdf", width = 16, height = 9)
VlnPlot(results_naive, "MA0101.1", pt.size = 0.1 )+ggtitle("RELA")+scale_fill_manual(values=BottleRocket3)
VlnPlot(results_naive, "MA1143.1", pt.size = 0.1 )+ggtitle("FOSL1::JUND")+scale_fill_manual(values=BottleRocket3)
VlnPlot(results_naive, "MA1512.1", pt.size = 0.1 )+ggtitle("KLF11")+scale_fill_manual(values=BottleRocket3)
dev.off()

## png 
##   2